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NUMERICAL METHODS AND COMPUTERS USED IN 
ELASTOH''DRODYNAMIC LUBRICATION 
Bernard J. Hamrock and John H. Tripp 

National Aeronautics and Space Administration 
Lewis Research Center 
Cleveland, Ohio 44135 

SUMMARY 

Some of the methods of obtaining approximate numerical solutions to 
boundary-value problems that arise in elastohyorodynamic lubrication are 
reviewed. The highlights of four general approaches (direct, inverse, qrasi- 
inverse, and Newton-Raphson) are sketched. Advantages and disadvantage's of 
these approaches are presented along with a flow chart showing some of the 
details of each. The basic question of numerical stability of the elasto- 
hydrodynamic lubrication solutions, especially in the pressure spike region, 
is considered. Computers used to solve this important class of lubrication 
problems are briefly described, with emphasis on recently developed 
supercomputers. 


INTRODUCTION 

Elastohydrodynamic lubrication is a form of fluid-film lubrication where 
elastic deformation of the lubricated surfaces becomes significant. It is 
usually associated with highly stressed machine components such as rolling- 
element bearings and gears. Historically, elastohydrodynamic lubrication may 
be viewed as one of the major developments in the field of lubrication in the 
twentieth century. Its recognition not only revealed the previously un- 
suspected regime of lubrication in highly stressed nonconformal machine ele- 
ments, but it also brought order to the complete spectrum of lubrication 
regimes, ranging from boundary to hydrodynamic. 

The present paper attempts to review the methods of obtaining approximate 
numerical solutions to boundary-value problems that a»"ise in tribology. The 
central task is to reduce the relevant differential and integral equations to 
algebraic ones that can be solved by familiar methods. Though the material 
is primarily illustrated by elastohydrodynamic lubrication problems, it is 
directly applicable to other areas of lubrication as well as to other engi- 
neering disciplines such a^ heat transfer and fluid mechanics. 

The highlignts of four main approaches to the elastohydrodynamic lubrica- 
tion problem, namely the direct method, the inverse method, the quasi-inverse 
method, and the Newton-Raphson method, are covered, and the advantages and dis- 
advantages of each method are discussed. The important question of numerical 
stability of solutions, especially in the pressure spike region, is considered. 
No attempt is made herein to be rigorous or complete. The chief purpose is to 
introduce the various techniques in a systematic manner and to indicate their 
general regimes of validity. Computers used to solve this important class of 
lubrication problem are briefly described with an emphasis on recently 
developed supercomputers. 
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Constants used to define density, tn /N 

O 

Influence coefficient used in elasticity calculation, m /N 
Modulus of elasticity, N/m^ 


Effective elastic modulus. 




Iteration damping factor 

Output/input amplitude ratio used in stability analysis 
Film shape, ti 

Film shape obtained from elasticity calculation, m 

Film shape obtained from inverse Reynolds equation, m 

Film shape constant, m 

Film thickness when dp/dx = 0, m 

Wavelength used in stability analysis 
2 

Pressure, N/m 

2 

Initial pressure, N/m*- 

2 

Inlet pressure, N/m 

2 

Pressure contact and outlet, N/m 
Effective radius, m 
Residual vector 
Curvature radius, m 
Banded matrices 

Mean velocity, (u^ + Uj^)/2, m/s 

Eigenvector 

Coordinate system 

2 

Pressure-viscosity coefficient of lubricant, m /N 

Elastic deformation, m 

Pressure difference, N/m^ 

Absolute viscosity at gage pressure, N s/m*^ 

2 

Viscosity at atmospheric pressure, N s/m 

Eigenvalue 

Poisson's ratio 

? n 

Lubricant density, N s /m"' 

2 4 

Density at atmospheric pressure, N s /m 
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Subscripts : 
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a Solid a 

b Solid b 

x,y Coordinate system 


BASIC EQUATIONS 

The elastohydrodynamic lubrication problem is as follows: given two 

elastic, independently rotating solids immersed in oil and pressed together by 
an external force, find the pressure distribution and film shape in the region 
of the lubricating contact. This requires calculating the pressure distribu- 
tion within the conjunction, at the same time allowing for the effects this 
pressure will have on the properties of the fluid and on the geometry of the 
elastic solids. The solution also provides the shape of the lubricant film, 
particularly the minimum film thickness between the solids. 

The basic equations used in elastohydrodynamic lubrication are as follows: 


Lubrication Equation (Reynolds Equation) 



where u = (u^ + %)I2 

Viscosity Variation 


n = 



( 2 ) 


where no ''s the coefficient of absolute dynamic viscosity at atmospheric 
pressure and a is the pressure-viscosity coefficient of the fluid. Values 
of ap may be as high as 10. 


Density variation 


Typical ly 


^ = *^0 1 +^bp) 

For mineral oils the values of a and b yield a maximum density increase of 
about 35 percent. 


6 = 


Elasticity Equation 
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h(x.y) = hg + ^ + 6(x,y) (5) 

where 

t 

1 ,J_.JL 

''x ''ax "bx 

1 =J_^J_ 

'^y '■ay ■'by 

Approaches that have been used in analyzing elastohydrodynamical ly lubricated 
conjunctions are presented briefly in the following sections. 

DIRECT METHOD 

The direct method solves the Reynolds equation for the pressure distribu- 
tion arising from a given film shape. A flow diagram of the direct method 
used by Hamrock and Dowson (1976) is presented in Fig. 1. The modified pres- 
sure shown on the left of Fig. 1 contains a damping factor to control numerical 
convergence. This pressure is then used to determine changes in the film 
shape. However, for maximum Hertzian contact stresses exceeding 0.5 GPa, even 
with this damping factor, the direct method is sometimes found to diverge. 
Nevertheless, Hamrock and Dowson (1981) were able to obtain useful formulas 
that cover a complete spectrum of contact geometries (ranging from point to 
line contacts), materials (hard and soft), and lubricant availability (fully 
flooded or starved conditions). These theoretical film thickness formulas 
were found to have a pleasing agreement with the experimental findings of 
Dalmaz and Godet (1978), Kunz and Winer (1977), and Koye and Winer (1980) even 
at relatively large maximum Hertzian contact stresses (1.5 GPa). The reason 
for this seems to originate partially from the linearity of the minimum film 
thickness - load relationship when plotted on a log-log scale. 

INVERSE METHOD 

The Reynolds equation (Eq. (1)) is normally regarded as determining the 
pressure distribution corresponding to a given film shape. In the inverse 
method adopted by Dowson and Higginson (1959), however, the equation is used 
to find the film shape responsible for the generation of a given pressure dis- 
tribution. That is, for a given load, an initial pressure profile (slightly 
different from Hertzian) is chosen. The film thickness is calculated twice, 
once by using the elasticity equation and again from the Reynolds equation. 

An elastohydrodynamic lubrication solution is obtained when the discrepancy 
between these two film profiles is sufficiently small. If the film shapes are 
not in agreement, the pressure profile is modified to improve the agreement 
between them. A flow diagram of the inverse method is shown in Fig. 2. 



Film Shape Equation 
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Some limitations of the inverse method are listed below: 

(1) It is not suited for lightly loaded cases where the film shape in 
the contact region is not parallel. 

(2) Although the computational method of Dowson and Higginson (1959) 
produced an acceptable solution in a small number of cycles, the procedure was 
not fully automatic. Judgement was needed in effecting the necessary modifi- 
cation of the pressure curve on the basis of the discrepancies between the 
elastic and inverse Reynolds film shape calculations. 

(3) It is only suitable for one-dimensional problems as pointed out by 
Rohde and Oh (1975). Where more than one spatial variable is involved, the 
Reynolds equation cannot be so readily integrated to express the film profile 
as a function of pressure. 


QUASI-INVERSF. METHOD 

With regard to limitation (3) given above, it should be pointed out that 
Evans and Snidle (1982) have recently used a quasi-inverse method in solving a 
two-dimensional elastohydrodynamic lubrication problem and have obtained good 
results for heavily loaded conditions. In contrast to solving an algebraic 
cubic equation as Dowson and Higginson (1959) did for the one-dimensional 
problem, Evans and Snidle (1982) have solved the Reynolds equation for the 
two-dimensional problem as a first-order uifferential equation in the film 
shape. They use a direct method in the inlet region, where the pressure in 
the conjunction is low, and the inverse method els'*where. Figure 3 shows the 
boundaries between the computing regions used in the quasi-inverse method. 

The circle in Fig. 3 represents the corresponding Hertzian area of dry contact. 
The curve AA is the starting point for the inverse solution of the Reynolds 
equation. The curve BB, which is downstream of AA, is the downstream boundary 
of the region over which a direct solution to the Reynolds equation is ob- 
tained. The procedure is summarized as follows: 

(1) Assume an initial pressure Po(x.y)* 

(2) Use elasticity to calculate film shape he(x,y). 

(3) Use the direct method to calculate pressure Pi(x,y) upstream of BB . 

(4) Adjust the film shape constant hQ until integrated load agrees with 
input load. 

Assume that tne pressure downstream of BB initially is Po(x,y). 

(5) Use the inverse hydrodynamic approach to calculate the film shape 
hi(x,y) downstream of AA. 

(6) Adjust pressure downstream of BB with 

r hg(x,y) 

P 2 (x,>) . Pp(x.y) 2 
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(7) Between AA and BB use a weighted mean of Pi(x,y) and P 2 (x,y) for 
the new pressure. 

(8) Continue calculations, returning to (2) until little change in pres- 
sure occurs throughout the conjunction. 

The position of BB in Fig. 3 is located where the viscosity-pressure exponent 
op « 1 or 2. The boundary AA is positioned two or more grid points upstream 
of BB. 

It appears that the Evans and Snidle (1982) approach is able to get con- 
verged solutions Hertzian contact stresses up to 1.5 GPa. This is a 
threefold increase in stress over the direct method and within the range that 
nonconformal contacts such as rolling element bearings and gears experience. 

A limitation of the approach is that a good initial pressure profile is needed 
since :he pressure in the zone between curves AA and BB is constrained to f 
very slowly by virtue of their very small separation. 

NEWTON-RAPHSON METHOD 

Both in its direct and inverse forms, the Reynolds equation is nonlinear 
(except the isoviscous, incompressible direct case), so that the Newton method 
along with its variations suggests itself as a powerful approach to the 
problem. Applied to an algebraic equation f(p) = 0, the Newton method pro- 
vides a systematic and usually convergent improvement on . ^ approximation 
p(*^^ to one of the roots. The simple algorithm that results can be expressed 
as p(*^‘*’l) = p('^) - f (p'’^0/f ' (p'^' ), which shows that the correction is always 
proportional to the residual error f{p(*^J). The Newton method for a func- 
tional equation likewise produces an algorithm for which the correction to the 
solution depends on the residual, with the chief difference being that the 
correction function itself satisfies a differential equation. With suitable 
boundary conditions this equation supplies convergent corrections to each 
previous approximation. 

Though the Newton-Raphson method is capable of solving the most general 
form of the Reynolds equation, we limit our demonstration of the method here 
to the one-dimensional, incompressible case. For this situation equation (1) 
can be written as 


f(p) = 


ax 


I t ax / ax 


The Newton-Raphson method implies that 




+ A 


df(p^^^) 

dp 


= 0 


where 


. . pi") 


A = p 


The index n refers to the n^*^ iterative approximation. But 


( 6 ) 


(7) 


( 8 ) 
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A 



+ 
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By making use of Eqs. (6) and (9), Eq. (7) can be rewritten as 


(9) 


0 



d 

37 


' . Ah^ dpt"> 
n dx 


^ ^ da 3h^6(a) dp^^^ 

n ^ n ax 

In Eq. (10) 6{a) Is the elastic deformation due to the pressure difference 

between the present and previous cycles at a given location. 


- 12u (10) 


Figure 4 is a flow diagram showing how the Newton-Raphson method applies 
to the elastohydrodynamic lubrication problem. There are three main loops. 
The innermost loop solves for the pressure difference a^*^ from Eq. (10) 
while the viscosity and film thickness at any x location are held constant. 
Within this loop the elastic deformation 5(a) due to the pressure difference 
A is continuously updated. With this converged solution on the pressure 
difference in the inner loop, the new pressure throughout the conjunction can 
be evaluated, and the viscosity and film shape are recalculated in the second 
loop. This loop is continued urtil there is little change in the pressure 
throughout the conjunction. The final loop requires that the load obtained 
from integrating the pressure is in agreement with the inputted load. 

Rohde and Oh (1975) applied higher order finite elements (cubic splines, 
cubic Hermite polynomials) with the Newton-Raphson iteration to solve lightly 
and moderately loaded cases (maximum Hertzian contact stress less than 
0.8 GPa). Numerical instabilities were found for larger maximum Hertzian 
contact stresses. 


CONVERGENCE AND STABILITY 

Hitherto several approaches to the elastohydrodynamic lubrication problem 
have been considered, each of wnich finds application in some regime of the 
operating parameters contained in the model equations. However, since none of 
the methods produces a solution (i.e., a compatible pressure distribution and 
film shape) in closed mathematical form, questions of convergence and stability 
of the various results assume considerable importance. Moreover, in arriving 
at an acceptable solution, it should not be overlooked that while uniqueness 
of solutions is generally assumed, it has never been rigorously demonstrated 
in the case of the particular nonlinearities presented by the elastohydro- 
dynamic lubrication problem. The rather weak and transparent basis for this 
assumption is simply that, as a matter of experience, only one numerical solu- 
tion is ever found for any specified set of operating parameters. 

Linear Analysis 

In view of the complications introduced by nonlinearity, combined with 
the great variety of numerical approximation schemes within the general cate- 
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gorles of finite element and finite difference, we do not attempt the most 
general analysis of convergence or stability here. Instead, we examine the 
behavior of a linearized method, such as the Newton-Raphson method, which can 
be suitably reduced to algebraic form to be solved by an Iterative scheme. 
Some comnents are also Included on the stability of the pressure distribution 
under conditions such that several extrema. Including a spike, may occur In 
the conjunction. 


Consider for example the innermost Iterative loop of the Newton method 
shown in the flow diagram of Fig. 4. The differential equation (Eq. (10)) for 
the pressure correction Is linear In both a and «, and as Implemented in 
Fig. 4 Its solution resembles that of time march methods In the transient 
problem with iteration count replacing time advance. An eigenvalue analysis 
of the convergence of such a loop is therefore appropriate. An equivalent 
linearized approach using an explicitly time-dependent perturbation was adopted 
by Kostreva (1983), who was able to construct a stability map In the space of 
two nondlmenslonal parameters. Curiously, however, most of the existing 
elastohydrodynamic solutions lie well within the unstable region of this map. 

Each term of Eq. (10) belongs to one of three types: (1) linear In a, 

(2) linear In a, or (3) independent of both a and 6. In the algebraic 
reduction of the derivatives invoked for types (1) and (2), their contribu- 
tions to Eq. (10) can be written in matrix form, respectively, as Tjjaj and 
Sij«j, where the summation convention is used. The matrices T and S will 
usually be banded with elements given by the initial (or current) approxima- 
tion to pressure p and film shape h. Single-subscripted quantities are 
vectors whose elements are the nodal values on the chosen discrete mesh. Type 

(3) terms are nonzero only when the pressure distribution fails to satisfy the 
Reynolds equation and hence are known ai. ^he vector of the residuals R-|. In 
this representation Eq. (10) becomes 


T, .A. + S. . - R. 
1j J iJ J 1 


If T is nonsingular, the direct solution of Eq. (11) Is 

A =* A + B6 


( 11 ) 


( 12 ) 


where A = T“^R and B = -T'^S. Since the elastic displacement vector 6 is 
not initially known, a is computed iteratively according to the prescription 




A + B«("l 


(13) 


Similarly for the inverse method the solution of Eq. (11) appears as follows: 


,(nn) 


« A' 


+ B'a 


(n) 


(13') 


where A' and B' require the nonsingularity of S. 

To complete the Iteration loop, the elastic equation relating a and 6 
Is needed. The elasticity equation such as Eq. (4) can be expressed as 

gin) . oa^"^ ( 14 ) 
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which Is not Iterative. The solution of Eq. (13), making use of Eq. (14). is 
thus 


. "f; (bo^a 

m>0 


(15) 


which shows explicitly the proportionality of a to the modified residual 
A s T"*R. The limiting value of a'”' becomes 


a - lim a^'’^ 


11m 

n4.» 


n-1 


m4) 


or 


a - (1 - BD)‘^A (16) 

The intermediate estimates a^^^ are just partial sums in the expansion of the 
inverse operator on the right side of Eq. (16). 

Similarly for the inverse method, making use of Eqs. (13') and (14) gives 

« - (1 - B'D"^)"^ A' (16‘) 

Because of the ill condition of 0, this se^ent of the inverse elastohydro- 
dynamic lubrication problem is never actually solved this way. Consequently 
the. inverse method discussed earlier ends up with two film shapes for each 
a^^^l instead of generating the usual pressure film shape loop. 

From Eq. (16) it is easy to see how the inner loop of Fig. 4 behaves. 

Let the eigenvector decomposition of the residual vector A be written as 

H V. (17) 


where A„ is the (scalar) coefficient of the eigenvector of the 

product matrix BD corresponding to eigenvalue Inserting Eq. (17) into 

Eq. (15) and summing the series yield 



(18) 


which converges if 1 x ^,1 <1 for all a. When this holds, ^qb. (16) and (18) 
both lead to 


a 


X )"^ X 
o' o 


(19) 
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This converged limit Is quickly attained If all iXa* *< 1* The most Impor- 
tant contributions to the sum are unfortunately the slowest to converge, and 
it may thus be useful to test for max ix^i . When this 1s close to 1, and 
particularly if it exceeds 1, the number of inner-loop cycles should be cur- 
tailed. Control is then passed to the second Iteration loop, where the matrix 
BD and the residual vector A are recomputed from the new pressure and its 
film shape. 


This behavior of the inner loop of the Newton-Raphson method remains 
essentially unaltered by the introduction of a damping factor, f. If the loop 
Is reentered with the weighted pressure distribution (1 - 
there is no change in the limiting iterate of Eq. (16). Only the intermediate 
steps by which this is reached are affected. It is worth noting, however, 
that any smallness criterion for the pressure changes used in determining an 
exit from the inner loop should, for consistency, be multiplied by (1 - f). 


Nonlinear System Response 


Useful as this analysis can be in carrying out the first iteration, it 
leaves unanswered the original question of overall convergence of the method. 
Specifically, it does not address the modification in the eigenvalues of BD 
when the middle iteration loop updates the pressure. Some limited appreciation 
of the stability at this. stage can be obtained by regarding this loop as an 
amplifier whose input p'^^ produces output p'*^ The parameters of the 
amplifier are dependent on the input, reflecting the nonlinearity or feedback 
aspects of elastohydrodynamic lubrication, and the objective is to achieve 
unit gain at any frequency by adjusting the parameters (matrix elements). 

Since the components of p are nodal pressures, frequency here refers to 
spatial oscillation in the solution vector. 

For illustrative purposes the treatment can be reduced to its simplest 
form by considering the same one-dimensional incompressible flow problem as in 
Eq. (6). The Reynolds equation in this situation can be integrated once to 
give 


, h - h 

^ ( 20 ) 

A typical solution of this equation is displayed in Fig. 5, which shows that 
over a length AB comparable to the Hertzian diameter the elastohydrodynamic 
lubrication film thickness may be substantially constant. The pressure and 
hence also the viscosity are large in this parallel film region, so that char- 
acteristics of the pressure profile are highly sensitive to the small differ- 
ence (h - hfn), where h^i is the film thickness at the pressure extrema. 

By superimposing a low-amplitude pressure ripple on the pressure input to 
the middle loop, the output ripple can be calculated with the aid of Eq. (20) 
and the one-dimension<J counterpart of Eq. (4). The output/ input amplitude 
ratio or gain g in a region of almost constant film thickness is then found 
to be given by the upproximate proportion 


10 


( 21 ) 


where K is the wave nutrtoer (reciprocal of wavelength) of the ripple. The 
constant of this approximate proportionality depends on the contac. geometry. 

Equation (21) shows qualitatively how the middle iteration loop is ex- 
pected to behave in the region AB of Fig. 5, where the pressure gradient is 
still small and where the pressure extrema may occur. For given u, no» and 
a values, h is essentially determined independently of load so that interest 
focuses on the p and K dependence of g. The larger g grows, the greater 
becomes the amplification of an)^ ripple, with consequent appearance of oscilla- 
tions in the pressure distribution. The equation confirms the general experi- 
ence that such are associated with high pressures, particularly on the upstream 
side of the spike in the vicinity of 6. It further suggests that oscillation 
bec(xnes of less concern at short wavelengths. Since the smallest wavelengths 
involved are of the order of the mesh -ize, numerical noise introduced by the 
algebraic representation of the basic equations or by finite ccxnputer accuracy 
is well attenuated. Greater difficulty is encountered at wavelengths compara- 
ble to the length of the parallel film region, and this ultimately may lead to 
distortion of much of the pressure curve there. 

Although the validity of Eq. (21) breaks down for still larger wave- 
lengths, it is clear that the system response to Such long-wavelength dis- 
turbance is given by its response to changes in ho, the additive constant in 
the film shape Eq. (5). This behavior is controlled separately in the third 
iteration loop of Fig. 4 and as such is not treated by present considerations. 

SUPERCOMPUTERS USED IN EHL 

The computers used to solve elaston.vdrodynamic lubrication problems in 
the 1950's were either mechanical or the first of the electronic computers and 
were considerably different from the supercomputers available today. The 
present-day powerful computing machines known as supercomputers have peak 
computing speeds exceeding 100 million operations per second. This is to be 
contrasted with the first commercial electronic computer delivered in 1951, 
which had a peak computing speed of around 600 operations per second. This 
means that the speed of large-scale scientific computers has doubled, on the 
average, every 2 years. Although the current performance levels of such 
machines owe much to the rapid advance of microelectronics, new concepts in 
computer architecture have been equally important. The term 'architecture' 
refers to the logical organization of the computer as it is seen by the pro- 
grammer. The architectural innovations of greatest significance are those 
that enable the machine to carry out many similar operations in parallel, or 
concurrently. This is referred to as vector computation. The latest super- 
computers allow the programmer to prescribe many different elementary steps to 
be executed simultaneously, whereas earlier computers could only handle a 
single sequence of elementary steps executed one at a time. This is referred 
to as scalar computation. Levine (1982) gives a thorough and useful descrip- 
tion of supercomputers. 

The most recent supercomputers provide a quantum jump in the speed of 
computation, and their effect in solving elastohydrodynamic lubrication prob- 
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1em$ In the next decade should be significant. Hamrock (i982) denx^nstrated 
the effect of using such con^uters In performing elastohydrodynamic film shape 
calculations. Both scalar and vector versions of the film shape code were 
presented. The film shape calculations were chosen since In elastohydrodynamic 
lubrication solutions over 80 percent of the computation t1n« Is spent In the 
straightforward calculation of elastic deformation. 

Before proceeding It might be well to point out more specifically what Is 
meant by scalar and vector coding. Scalar coding Is conventional sequential 
coding that has been ured universally for digital computers since their 
earliest development. the other hand, the process of writing programs In 
vector code consists of organizing programs so that the maximum possible num- 
ber of arithmetic and logical operations can be processed In parallel. Such 
'vectorizatlon' Is most easily Introduced and Illustrated by an example. 

Consider arrays A and B, each consisting of 100 nuntf)ers, which are to 
be added to form array C, where Cj •» Aj + Bj, 1 ■ I, . . 100. The tradi- 

tional 'scalar* computer executes five assembly language Instructions each 
100 times. There are two reads from memory (Ai and B-j), one addition, one 
store to memory {for C^), and an Instruction that increments a counter, tests i 
and branches back to load the next pair of input operands. Thus 500 scalar 
Instructions are executed In consecutive sequence to add arrays A and b. 

A vector computer's compile*" can generate 'vector' object code, which executes 
very differently. The vector code for adding the 100 pairs of operands still 
uses the first four of these In;. tractions for each pair. But instead of con- 
tinually branching back for the next pair, execution proceeds by continuously 
streaming operands from central memory Into the central processor, where the 
addition takes place, and continuously streaming answers back Into fromory In 
'pipeline' fashion. During execution of the vectorized addition son« el^nents 
of A and B are being read from central memory, some elements of A and 
B are undergoing addition in the vector pipeline, and some are being written 
in memory. 

This brief example illustrates the major conceptual difference between 
the two modes of operation. It should . clear from this that many kinds of 
repeated instructions can be processed in vector form, with consequent enormous 
savings in overall time. It is also apparent that there are other situations 
in which the pipeline mode will not operate. Perhaps the most common of these 
in elastohydrodynamic Rubrication problems arise when a recursive or iterative 
process requires the i result before it can generate the (i + 1)^", as in 
direct integration of a differential equation. Often a way may be found to 
avoid such conditions, but even when this is not possible and trie computer is 
forced into the scalar mode, its operation speed is still considerably superior 
to previous generation computers. 

Table I shows the scalar and vector computation times expressed in nano- 
seconds for varying vector 'length', which here corresponds to the number of 
nodes In the finite difference representation of the elastohydrodynamic lubri- 
cation solutions. A considerable decrease in computation time is seen for the 
vector confutation over that of the scalar computation. The last column of 
the table gives the ratio of scalar to vector times and indicates that, as the 
vector ength increases, the ratio also increases significantly, although 
eventually an asymptotic condition is approached. The run times for the two 
types of code indicate that over a 50-to-l speedup of scalar to vector compu- 
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tatlon time is obtained for vector lengths typically used in elastohydrodynamic 
lubrication analysis. 


CONCLUSIONS AND OUTLOOK 

The aim herein has been to examine the principal approaches devised to 
tackle the technologically Important elastohydrodynamic lubrication problem. 
Without entering into details of numerical approximation schemes applicable to 
the nonlinear equations involved, it is apparent that under the general head- 
ings of finite difference and finite element, sufficiently powerful methods 
exist to produce some useful solutions to the problem. Executing such schemes 
by a multigrid approach could lead to further substantial improvements over 
the conventional fixed-mesh methods most widely used at present. The avail- 
ability of supercomputers combining high speed and large memory capacity 
removes most obstacles to implementing the chosen scheme, so that the effi- 
ciency or convergence rate of the method is no longer the overriding concern. 
The question of central importance has now become the stability of the solu- 
tion method. 

Experience with finite element techniques in elastohydrodynamic lubrica- 
tion applications is accumulating, as evidenced by the two sessions of this 
conference devoted to them, and it is likely that such techniques could have 
some advantages over the currently more familiar finite difference methods. 
Thus, it often appears that significantly fewer nodal values are required in 
the finite element representation than for finite differences. Nevertheless 
the strong impetus behind finite element analysis, namely, its suitability for 
problems of little or no symmetry, as, for example, in classical elastostat Ics, 
is lacking in elastohydrodynamic lubrication, where nominal bearing geometry 
is often required to be highly syirenetrical . 

Practical experience indicates an advantage to using an inverse method at 
higher maxiriium pressures and conversely a direct method at lower pressures. 

The linearized Newton-Raphson method of improving solutions to functional 
equations runs in either the direct or inverse mode according to which ap- 
proach to the hydrodynamic equation is taken. So far, no serious attempt has 
been made to invert the elasticity equation, but Saint-Venant's principle 
strongly suggests that such an approach would be hopeless. 

Turning to questions of stability, again the issue is not between finite 
difference and finite element. Both approaches can be analyzed in the same 
manner. A beginning at such analysis has been made by Kostreva (1983), whose 
stability map uses a criterion based on the eigenvalue behavior of a certain 
linearized form of the Reynolds equation. 

A feature shared by many elastohydrodynamic lubrication calculations has 
been that their implementation by computer is not fully automatic. The solu- 
tion evolves stepwise from the initial guess, and the criterion for passing to 
the next step (iteration) requires some judgement or a further guess not easily 
reduced to a sequence of simple branching decisions. As such the pathway to a 
final solution has been neither completely described nor properly optimized 
and the nonquantifiable skills of the programmer have played a prominent role. 
To exploit more fully the advantages of automatic computation, it is thus 
necessary to develop improved understanding of the convergence properties of 
the various steps and the ability to decide if the current step is actually 
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moving toward or away from a solution. Since several constraints are imposed 
which it is difficult to tighten simultaneously, the process at intermediate 
stages may seem to be diverging from the viewpoint of some single constraint. 
It is common, for example, to find a pressure profile making improvements in 
the film shape while diverging from the required load, to which adjustments 
are made only at a later stage. An optiml 2 ed strategy might handle the load 
at the same time. In brief, criteria are needed to enable control to pass 
automatically between the iterative loops. 

It was to this end that the analysis of the behavior of one type of 
iteration loop was undertaken here. A criterion for convergence in terms of 
eigenvalues of the operator representing the net effect of the loop was 
derived. Although it still remains to develop this into concise executable 
form, progress toward this goal is being made and further effort in this 
direction is certainly worthwhile. 

Similarly, the treatment given here of the much more difficult stability 
problem associatfd wHh the high-pressure region of the contact is a beginning 
attempt to understand the conditions for which reliable solutions can be gen- 
erated there. Further work, however, is needed to improve the analysis as the 
spike region itself is entered, where the pressure gradient grows large. 

If the potential of such approaches to program automation is fulfilled, 
then prospects are excellent for considering additional aspects of the elasto- 
hydrodynamic lubrication problem, such as 

(1) Non-Newtonian effects of the fluid 

(21 Surface roughness effects 

(3) Thermal effects 

With the use of supercomputers for the incorporation of these coupled effects 
it can be predicted that run times will be dramatically reduced during this 
decade. Such invest! "tions will then become much more feasible than in the 
past, with a correspon <ng gain in the reliability and range of application of 
the results. 
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TABLE I. - COMPARISON OF SCALAR AND VECTOR 
COMPUTER RUN TIMES FOR DIFFERENT 
VECTOR LENGTHS 


liill 

Scalar time, 
nsec 

Vector time, 
nsec 

Scalar time 
Vector time 

10 

241 

70 

3.44 

30 

1 179 

155 

7.61 

50 

2 726 

243 

11.22 

100 

9 252 

483 

19.28 

300 

73 357 

1 933 

37.95 

500 

198 282 

4 188 

47.35 

700 

383 997 

7 242 

53.02 

1 000 

776 588 

13 363 

58.11 

3 000 

6 890 354 

99 988 

68.91 

5 000 

19 084 530 

266 602 

71.58 

7 000 

37 358 446 

513 218 

72.79 

10 000 

76 171 860 

1 033 169 

73.73 
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Figure 3, - Contact area showing boundaries between computing 
regions used in quasi -inverse method. (From Evans and 
Snidle U982)) 
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Figure 4 - Flow diagram of Newton-Raphson method. 
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Figure 3. - example of an ehl pressure distribution and film profile. The 
operating parameters have been chosen to produce a minimum in the 
pressure curve between the Hertzian maximum and the spike near the 
exit Between points A and B, the film thickness is nearly constant - 
as a consequence of the similarity of the pressure to the Hertzian dis- 
tribution. shown for comparison. 



